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A limited mobility nonequilibrium solid-on-solid dynamical model for kinetic surface growth is in- 
troduced as a simple description for the morphological evolution of a growing interface under random 
vapor deposition and surface diffusion bias conditions. Simulations using a local coordination depen- 
dent instantaneous relaxation of the deposited atoms produce complex surface mound morphologies 
whose dynamical evolution is inconsistent with all the proposed continuum surface growth equations. 
For any finite bias, mound coarsening is found to be only an initial transient which vanishes asymp- 
totically, with the asymptotic growth exponent being 0.5 in both 1+1 and 2+1 dimensions. Possible 
experimental implications of the proposed limited mobility nonequilibrium model for real interface 
growth under a surface diffusion bias are critically discussed. 
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An atom moving on a free surface is known to en- 
counter an additional potential barrier, often called a 
surface diffusion bias Q , as it approaches a step from the 
upper terrace — there is no such extra barrier for an atom 
approaching the step from the lower terrace (the surface 
step separates the upper and the lower terrace). Since 
this diffusion bias makes it preferentially more likely for 
an atom to attach itself to the upper terrace than the 
lower one, it leads to mound (or pyramid) - type struc- 
tures on the surface under growth conditions as deposited 
atoms are probabilistically less able to come down from 
upper to lower terraces. This dynamical growth behavior 
is sometimes called an "instability" because a flat ( "sin- 
gular" ) two dimensional surface growing under a sur- 
face diffusion bias is unstable toward three dimensional 
mound/pyramid formation. There has been a great deal 
of recent interest p)-|2l|] in the morphological evolution 
of growing interfaces under nonequilibrium growth con- 
ditions in the presence of such a surface diffusion bias. 
In this paper we propose a minimal nonequilibrium cel- 
lular automata - type atomistic growth model for ideal 
molecular beam epitaxial - type random vapor deposition 
growth under a surface diffusion bias. Extensive stochas- 
tic simulation results presented in this paper establish 
the morphological evolution of a surface growing under 
diffusion bias conditions to be surprisingly complex even 
for this extremely simple minimal model. Various criti- 
cal growth exponents, which asymptotically describe the 
large-scale dynamical evolution of the growing surface in 
our minimal discrete growth model, are inconsistent with 
all the proposed continuum theories for nonequilibrium 
surface growth under diffusion bias conditions. Our re- 
sults based on our extensive study of this minimal model 
lead to the conclusion that a continuum description for 
nonequilibrium growth under a surface diffusion bias does 
not exist (even for this extremely simple minimal model) 
and may require a theoretical formulation which is sub- 
stantially different from the ones currently existing in 
the literature. Our results in the initial non- asymptotic 
transient growth regime (lasting upto several hundred or 



a few thousand layers of growth) do, however, agree with 
existing theoretical and (many, but not all) experimental 
findings in the literature. 

In Fig. 1(a) we schematically show our solid-on-solid 
(SOS) nonequilibrium growth model : (1) Atoms are de- 
posited randomly (with an average rate of 1 layer/unit 
time, which defines the unit of time in the growth prob- 
lem — the length unit is the lattice spacing taken to 
be the same along the substrate plane and the growth 
direction) and sequentially on the surface starting with 
a flat substrate; (2) a deposited atom is incorporated 
instantaneously if it has at least one lateral nearest- 
neighbor atom; (3) singly coordinated deposited atoms 
(i.e. the ones without any lateral neighbors) could in- 
stantaneously relax to a neighboring site within a dif- 
fusion length of / provided the neighboring site of in- 
corporation has a higher coordination than the orig- 
inal deposition site; (4) the instantaneous relaxation 
process is constrained by two probabilities Pt, and Ptj 
(0 < Pl,Pu < 1) where Pl(u) 1S the probability for the 
atom to attach itself to the lower(upper) terrace after 
relaxation (note that a "terrace" here could be just one 
other atom). The surface diffusion bias is implemented 
in our model by taking Pjj > Pt,, making it more likely 
for atoms to attach to the upper terrace. Under the sur- 
face diffusion bias, therefore, an atom deposited at the 
top of a step edge feels a barrier (whose strength is con- 
trolled by Pu/Pl) in coming down compared with an 
atom at the lower terrace attaching itself to the step. 
Our model is well-defined for any value of the diffusion 
length I including the most commonly studied situation 
of nearest-neighbor relaxation (I = 1). (We should em- 
phasize, however, that the definition of a surface diffusion 
bias is not unique even within our extremely simple lim- 
ited mobility nonequilibrium growth model — what we 
study in this Letter is the so-called edge diffusion bias 
p+[.) We have carried out extensive simulations both in 
1 + 1 and 2 + 1 dimensions (d) varying Pt, , Pjj as well as 
I, also including in our simulations the inverse situation 
(the so-called 'negative' bias condition) with Pt, > Pjj so 
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FIG. 1. (a) The mound evolution in 1+1 dimensions show- 
ing a section of 500 middle lattice sites from a substrate size 
of L = 10000 at 10 - 10 6 ML (Pu = 1; Pl = 0.5). Lower 
inset: Schematic configuration defining growth rules in 1+1 
dimensions, (b) The details of the initial mound evolution in 
(a) for 10— 10 4 ML, showing that the mounds form very early 
in the growth. Somewhat noisy mounds can already be seen 
at 10 ML, and the coarsening of the mounds between 10— 10 4 
ML is clearly seen in the results. (c,d) The d=2+l growth 
morphologies on a 100 x 100 substrate [Pu = 1; Pl ~ 0.5) at 
(c) 10 3 ML, and (d) 10 6 ML. 



that deposited atoms preferentially come down attach- 
ing themselves to lower steps producing in the process a 
smooth growth morphology. Because of lack of space we 
do not present here our 'negative' bias (1 > Pl > Pu) re- 
sults (to be published elsewhere) except to note that the 
smooth dynamical growth morphology under our nega- 
tive bias model obeys exactly the expected Q linear Ed- 
wards - Wilkinson universality |pl| . Our growth model 
is the most obvious finite bias generalization of the well- 
studied cellular automaton model referred to as the DT 
model or the 1+ model in the literature pl|-p3[|. 

Before presenting our numerical results we point out 
two important features of our growth model : (1) For 
Pl = Pu = 1 our model reduces to the one introduced 
in ref. 22, often called the DT model, (and studied ex- 
tensively pl|-p6| in the literature) as a minimal model 
for molecular beam epitaxy in the absence of any diffu- 
sion bias ; (2) w e find, in complete agreement with earlier 
findings |22],[23|] in the absence of diffusion bias, that the 
diffusion length I is an irrelevant variable (even in the 
presence of bias) which does not affect any of our calcu- 
lated critical exponents (but does affect finite size correc- 
tions — increasing I requires a concomitant increase in 
the system size to reduce finite size effects). To demon- 
strate this, we compare two systems (see Fig. 2(c)), one 
with diffusion length I = 1 while the other with 1 = 5. 
Both systems are on sufficiently large substrates of size 
L = 10 4 to prevent finite size effect and both have the 
same bias strength (Pu = 1.0 and Pl = 0.9). Except for 
the expected layer-by-layer growth seen in the first few 
layers in the I = 5 system, we see the same critical behav- 
ior from the two systems with the coarsening processes 
slightly faster in the I = 5 system in early time regime. 
In the rest of this paper (except for Fig. 2(c) where we 
present some representative I — 5 results), we present 
our I = 1 simulation results emphasizing that our crit- 
ical exponents are independent of I provided finite size 
effects are appropriately accounted for. Our calculated 
exponents are also independent of the precise values of 
Pu and P L ( < Pu ) as found |2|j2^] in the unbiased 
Pu = Pl case. 

In Fig.l we show our representative d=l+l (a and 
b) and 2+1 (c and d) simulated dynamical growth mor- 
phology evolution. The diffusion bias produces mounded 
structures which are visually statistically scale invari- 
ant only on length scales much larger (or smaller) than 
the typical mound size. Note that the mounding in the 
growth morphology starts very early during growth and 
is already prominent in the first 100 monolayers (ML) as 
is obvious in Fig. 1(b). In producing our final results 
we utilize a noise reduction |2q] technique which accepts 
only a fraction of the attempted kinetic events, and in 
the process produces smoother results (reducing noise ef- 
fects) without affecting the critical exponents. We have 
explicitly checked that the noise reduction technique does 
not change our calculated critical exponents |^8| and has 
only the cosmetic effect of suppressing noise in our sim- 
ulated growth morphology. 
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To proceed quantitatively we now introduce the dy- 
namic scaling ansatz 21 - 23 which seems to describe well 
all our simulated results. We have studied the root mean 
square surface width or surface roughness (W), the av- 
erage mound size (R), the average mound height (H), 
and the average mound slope (M) as functions of growth 
time. We have also studied the various moments of dy- 
namical height-height correlation function, and these cor- 
relation function results (to be reported elsewhere) are 
consistent with the ones obtained from our study of W(t), 
R(t), H(t), and M(t). The dynamical scaling ansatz in 
the context of the evolving mound morphologies can be 
written as power laws in growth time (which is equivalent 
to power laws in the average film thickness) : W(t) ~ t@ 
; R(t) ~ t n ; H(t) ~ t K ; M(t) ~ t x ; £(i) - t^ z , where 
is the lateral correlation length (with z as the dy- 
namical exponent) and /3, n, K, A, z are various growth 
exponents which are not necessarily independent. We 
find in all our simulations n ~ z , and thus the coars- 
ening exponent n, which describes how the individual 
mound sizes increase in time, is the same as the inverse 
dynamical exponent in our model. We also find = k 
in all our results, which is understandable in a mound- 
dominated morphology. In addition, all our results sat- 
isfy the expected exponent identity (3 = n = n+X because 
the mound slope M ~ H/R. The evolving growth mor- 
phology is thus completely defined by two independent 
critical exponents f3 (the growth exponent) and n (the 
coarsening exponent), which is similar to the standard 
(i.e. without any diffusion bias) dynamic scaling situa- 
tion |2l|-|23j] where and z (= n^ 1 in the presence of 
diffusion bias) completely define the scaling properties. 

In Figs. 2 and 3 we show our representative scaling 
results in d=l+l (Fig. 2) and 2+1 (Fig. 3) for nonequi- 
librium growth under surface diffusion bias conditions. 
It is clear that we consistently find the growth expo- 
nent (3 ~ 0.5 in both d=l+l and 2+1 in the long time 
asymptotic limit independent of Pl and Pjj as long as 
Pl/Pjj < 1- This j3 ~ 0.5 is, however, different from 
the usual Poisson growth under pure random deposition 
with no relaxation where there are no lateral correla- 
tions. Our calculated asymptotic coarsening exponent n 
in both d=l+l and 2+1 is essentially zero (< 0.1) at long 
times. In all our results we find the effective coarsening 
(growth) exponent showing a crossover from n (/?) « 0.2 
(0.25 or 0.33 depending on d=2+l or 1+1) at early times 
(1 < t < 10 3 ) to a rather small (large) value ( n < 0.1, 
j3 k 0.5 as i oo) at long times — we believe the 
asymptotic n (fj) to be zero (half) in our model. Our 
calculated steepening exponent A satisfies the exponent 
identity A = /3 — n rather well, indicating that steepening 
(M ~ £ A ) and coarsening (R ~ t n ) are competing pro- 
cesses. We find that during the initial transient regime 
(for 1 < t < 10 3 depending on d and Pl/Pjj), when con- 
siderable mound coarsening takes place (n ~ 0.2), A is 
rather small (A < 0.1) and does not change much. We 
note that the mound formation dominates our growth 
morphology even during the initial transient — mound- 
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FIG. 2. (a) The surface roughness W in the L = 10000 
system with Pu = 1; Pl ~ 0.5. Left inset: the growth 
exponent f3 calculated from the local derivative of log 10 W 
with respect to log 10 t. Right inset: Average mound size vs 
time, (b) W in L = 10000 systems with Pu = 1;P L = 0.25 
(upper curve) and Pu = 1; Pl ~ 0.75 (lower curve). Left 
inset: Average mound size vs time for the system with 
Pu = 1; Pl = 0.25. Right inset: The same plot for 
Pu = 1; Pl = 0.75. (c) Results for diffusion length / = 1 
compare with I — 5 (all other results in this paper correspond 
to I = 1). The results for I — 1 and I — 5 are parallel to 
each other for both W and R (inset) as a function of time, 
indicating exponent universality with respect to the diffusion 
length. In this figure, Pu = 1.0; Pl = 0.9. The bending in 
the W(t), R(t) curves in log — log plots shown here indicates 
the crossover behavior from the initial transient regime to the 
asymptotic regime as discussed in the text and shown in Fig. 
2 (a) and (b). The layer- by- layer oscillation can clearly be 
seen upto 10 (3) ML for 1 = 5(1 = 1). 
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FIG. 3. (a) The surface roughness W in the 100 x 100 
system with Pu = 1; Pj, = 0.5. Left inset: the local growth 
exponent (3. Right inset: Average mound height vs time, (b) 
The average mound size in the same 100 x 100 system. Left 
inset: the local coarsening exponent n calculated from the 
local derivative of log 10 R with respect to log 10 t. Right inset: 
the average mound slope vs time. 

ing starts early, coarsens rapidly, and then coarsening 
slows down or almost stops. After the initial transient, 
however, A is finite and large (A > 0.4), indicating sig- 
nificant steepening of the mounds. Based on our sim- 
ulation results we are compelled to conclude that the 
evolutionary behavior of our mound dynamics makes a 
crossover (at t = t* < 10 3 ) from a coarsening-dominated 
preasymptotic regime / (nj ~ 0.2, A/ < 0.1) to a 
steepening-dominated regime II [nn < 0.1, \n > 0.4) 
which we believe to be the asymptotic regime. The ini- 
tial transient (regime /: t < t*) can be construed as a 
"slope selection" regime where the steepening exponent 
A is very small (and the coarsening exponent, n ~ 0.2, 
is approximately a constant), but the asymptotic long 
time behavior (regime //: t > t*) is clearly dominated 
by slope steepening (large A) with coarsening essentially 
dying down (n < 0.1). The crossover time t* between 
regime / (coarsening) and // (steepening) depends on 
the details of the model (e.g. d, L, Pl, Pjj, I), and could 
be quite large (t* ~ 10 2 - 10 4 ). 

We emphasize that our initial transient (t < t*) expo- 
nents in regime I ((3j ~ 0.26; nj ~ 0.17; A/ ~ 0.09 in 
d=2+l, Fig. 3) agree quantitatively with several [p|,p|p5| 
non-minimal detailed (temperature-dependent Arrhcnius 



diffusion) growth simulations as well as with direct nu- 
merical simulations jl(i(| of a proposed (empirical) con- 
tinuum growth equation, which have all been claimed 
^|,|6| Jl^Jl5[] to be in good agreement with the observed 
coarsening behavior in various epitaxial growth exper- 
iments. Other simulations and theories PUT?]] find no 
slope selection, which also agree with some experiments. 
We believe that the observed slope selection (A rj 0) in 
experiments and simulations is a (long-lasting) transient 
(our regime I, t < t*) behavior, which should disappear 
in the asymptotic regime, at least in models without any 
artificial transient downward mobility. Theoretically, a 
distinction has been made p6|Jl7|JTs|^H between 'weak' 
and 'strong' diffusion bias cases, corresponding respec- 
tively to our regime / ('mounding') and II ('steepen- 
ing'), respectively. Our results show that this theoreti- 
cal distinction is not meaningful because the 'weak' bias 
case crosses over to the 'strong' bias case for t > t*, 
and the long time asymptotic regime is invariably the 
'strong' bias steepening regime for any finite diffusion 
bias. (There could be accidental slope selections at rather 
large slopes when crystallographic orientations are taken 
into account |l^,|2l|, a process neglected in our minimal 
growth model.) 

In comparing with the existing continuum growth equa- 
tion results we find that none can quantitatively explain 
all our findings. Golubovic |Q predicts f3 = 0.5, which 
is consistent with our asymptotic result (/3/j ~ 0.5), but 
his finding of n = A = 1/4 in both d=l+l and 2+1 is 
inconsistent with our asymptotic results (asymptotically 
n < 0.1, and A ~ 0.4 — 0.5) while being approximately 
consistent with our initial transient results (for t < t*). 
The analytic results of Rost and Krug also cannot 
explain our results, because they predict, in agreement 
with Golubovic, that if /3 = 1/2, then n = A = 1/4. 
We also find our asymptotic (3 to be essentially 0.5 in- 
dependent of the actual value of Pu / Pl , which disagrees 
with ref. 19. Interestingly our initial transient regime 
is, in fact, approximately consistent with the theory of 
ref. 19. We have approximate slope selection (A < 0.1) 
only in the initial transient regime, which could, however, 
be of considerable experimental relevance because the 
pre-asymptotic regime is a long-lasting transient. The 
only prior work in the literature that has some similarity 
to our results is that of Villain and collaborators PS, 
who in a one dimensional deterministic (i.e. without the 
deposition beam shot noise) macroscopic continuum de- 
scription of growth roughness and instabilities under sur- 
face diffusion bias, called the Zeno model by the 
authors, found, in agreement with our atomistic two di- 
mensional stochastic cellular automaton simulation re- 
sults, a scenario in which coarsening becomes "extremely 
slow after the mounds have reached a (characteristic) 
radius" which is reminiscent of our crossover from 
weak bias (i < t*) to strong bias (t > t*). While the 
precise relationship between our two (and one) dimen- 
sional atomistic/stochastic model and their one di- 
mensional macroscopic/deterministic model is unclear at 
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the present time, it is interesting to note that the au- 
thors of refs. came to a similar negative conclusion as 
we do about the non-existence of any continuum growth 
equation describing the strong bias asymptotic regime 
where the microscopic lattice size may play a crucial role 
p5fl . Finally, we note that a very recent experimental 
work |2^] reports growth morphology of Ge(OOl) surface 
which agrees qualitatively with the scenario predicted in 
this paper, namely, that even for a very weak diffusion 
bias, the mound slope continues to increase without any 
observable coarsening. 

Although our diffusion bias model is an extremely sim- 
ple limited mobility model which may be viewed as unre- 
alistic, our d=l+l results are remarkably similar to those 
obtained from a study of a full temperature-dependent 
Arrhenius hopping model with step edge barrier 
This study jll| offered essentially the same picture as 
what we present here, i.e. a smaller value of the growth 
exponent that crosses over to (3 ss 0.5 at larger time, 
and a very large dynamical exponent corresponding to 
very little coarsening (n — > 0) after approximately 100 
monolayers of deposition. To be specific, the dynamical 
exponent determined by the growth of the correlation 
length is z f=a 16.6 in the temperature-dependent edge 
bias model jll| while our study yields z > 10. Although 
the dynamical exponents from the two studies are not 
exactly the same, they are both exceptionally large in- 
dicating the mound coarsening process to be negligible 
in the large time regime. This qualitative agreement be- 
tween our simple limited mobility results and a d=l+l 
full diffusion results argue strongly in favor of our mini- 
mal growth model being of reasonably general qualitative 
validity in experimental situations. 

Finally, we note that the introduction of limited mo- 
bility models ]22| ] has opened a whole new way of study- 
ing kinetic surface roughening in molecular beam epi- 
taxial growth. These limited mobility nonequilibrium 
models make it possible to study very large systems in 
very large time limit when it is impossible to do so in 
the realistic temperature-dependent full- activated diffu- 
sion models. Particularly, we point out the success of 
the DT model |^2| in providing an excellent zeroth order 
description of molecular beam epitaxial growth in the ab- 
sence of any surface diffusion bias. The d=2+l critical 
exponents (2(| in the unbiased (DT) model p3|, which 
belongs to the same universality class as the conserved 
fourth-order nonlinear continuum MBE growth equation 

HHl, are f 3 = °' 25 ~ °' 2 and " - 0-6 - 0.7, which are 
in quantitative agreement with a number of experimental 
measurements [|2l) where surface diffusion bias is thought 
to be dynamically unimportant. With this in mind, it 
is, therefore, conceivable that our study of the limited 
mobility model which includes surface diffusion bias (i.e. 
a generalized version of the DT 22 model) presented 
in this paper may benefit the subject in a way similar 
to what the DT model did for the unbiased molecular 
beam epitaxy growth study. This is particularly signifi- 
cant since there are still many open questions regarding 



the interface growth under surface diffusion bias condi- 
tion. Since the model we study here is a generalized DT 
model 1 22 1 and an approximate continuum description for 
the original unbiased model [^2| has recently been devel- 
oped pq] , one could use that as the starting point to con- 
struct a continuum growth model for the biased growth 
situation. Such a continuum description is, however, ex- 
tremely complex p5| as it requires the existence of an in- 
finite number of nonlinear terms in the growth equation, 
and it therefore remains unclear whether a meaningful 
continuum description for our discrete simulation results 
is indeed possible @,|||]]. 

In conclusion, we want to emphasize the fact that the 
limited mobility model studied in this paper is an ex- 
tremely simple model ("the minimal modeT in the sense 
that it is perhaps the simplest nonequilibrium model 
which captures the minimal features of growth under sur- 
face diffusion bias), and realistic growth under experi- 
mental conditions should be substantially more complex 
than this minimal model (we speculate that this mini- 
mal model is in the same growth universality class as 
realistic growth under a surface diffusion bias for reasons 
discussed above, but we certainly cannot prove that at 
this incomplete stage of the development of the subject. 
A word of caution is in order in comparing experimental 
results with our calculated critical exponents because of 
the extreme simplicity and the limited mobility nature of 
our nonequilibrium growth model.) If we do not have a 
reasonable theoretical understanding (from a continuum 
equation approach) of even such a simple minimal model, 
as we have found in this paper, then current efforts |l] |2lJ 
at understanding realistic growth under surface diffusion 
bias must be quite futile. The main weakness of limited 
mobility model (of the type presented here) is that they 
are manifestly nonequilibrium model — this, however, 
should not be a particularly serious problem in the con- 
text of the mound/pyramid formation in the growth mor- 
phology, which by definition is a nonequilibrium effect 
and must disappear in a properly equilibrated surface. 
Our conclusion based on the results presented in this 
Letter is that a continuum growth equation for nonequi- 
librium growth under a surface diffusion bias does not 
exist at the present time. 

This work is supported by the US-ONR and the NSF- 
DMR-MRSEC. 
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